Проект 2
Подготовка к работе
Нам понадобятся следующие пакеты.
2. Подготовим данные
Проверим отсутсвующие значения и выбросим их при необходимости.
numberOfNA <- length(which(is.na(Boston)==T))
if(numberOfNA>0) {
Boston <- Boston[complete.cases(Boston),]
}Создадим данные для тестирования модели.
set.seed(123)
split <- sample.split(Boston,SplitRatio = 0.75) #assigns booleans to a new coloumn based on the split ratio
train <- subset(Boston,split==TRUE)
test <- subset(Boston,split==FALSE)3. Разведочный анализ.
Посмотрим структуру датасета и базовую статистику.
str(Boston)## 'data.frame': 506 obs. of 14 variables:
## $ crim : num 0.00632 0.02731 0.02729 0.03237 0.06905 ...
## $ zn : num 18 0 0 0 0 0 12.5 12.5 12.5 12.5 ...
## $ indus : num 2.31 7.07 7.07 2.18 2.18 2.18 7.87 7.87 7.87 7.87 ...
## $ chas : int 0 0 0 0 0 0 0 0 0 0 ...
## $ nox : num 0.538 0.469 0.469 0.458 0.458 0.458 0.524 0.524 0.524 0.524 ...
## $ rm : num 6.58 6.42 7.18 7 7.15 ...
## $ age : num 65.2 78.9 61.1 45.8 54.2 58.7 66.6 96.1 100 85.9 ...
## $ dis : num 4.09 4.97 4.97 6.06 6.06 ...
## $ rad : int 1 2 2 3 3 3 5 5 5 5 ...
## $ tax : num 296 242 242 222 222 222 311 311 311 311 ...
## $ ptratio: num 15.3 17.8 17.8 18.7 18.7 18.7 15.2 15.2 15.2 15.2 ...
## $ black : num 397 397 393 395 397 ...
## $ lstat : num 4.98 9.14 4.03 2.94 5.33 ...
## $ medv : num 24 21.6 34.7 33.4 36.2 28.7 22.9 27.1 16.5 18.9 ...summary(Boston)## crim zn indus chas
## Min. : 0.00632 Min. : 0.00 Min. : 0.46 Min. :0.00000
## 1st Qu.: 0.08204 1st Qu.: 0.00 1st Qu.: 5.19 1st Qu.:0.00000
## Median : 0.25651 Median : 0.00 Median : 9.69 Median :0.00000
## Mean : 3.61352 Mean : 11.36 Mean :11.14 Mean :0.06917
## 3rd Qu.: 3.67708 3rd Qu.: 12.50 3rd Qu.:18.10 3rd Qu.:0.00000
## Max. :88.97620 Max. :100.00 Max. :27.74 Max. :1.00000
## nox rm age dis
## Min. :0.3850 Min. :3.561 Min. : 2.90 Min. : 1.130
## 1st Qu.:0.4490 1st Qu.:5.886 1st Qu.: 45.02 1st Qu.: 2.100
## Median :0.5380 Median :6.208 Median : 77.50 Median : 3.207
## Mean :0.5547 Mean :6.285 Mean : 68.57 Mean : 3.795
## 3rd Qu.:0.6240 3rd Qu.:6.623 3rd Qu.: 94.08 3rd Qu.: 5.188
## Max. :0.8710 Max. :8.780 Max. :100.00 Max. :12.127
## rad tax ptratio black
## Min. : 1.000 Min. :187.0 Min. :12.60 Min. : 0.32
## 1st Qu.: 4.000 1st Qu.:279.0 1st Qu.:17.40 1st Qu.:375.38
## Median : 5.000 Median :330.0 Median :19.05 Median :391.44
## Mean : 9.549 Mean :408.2 Mean :18.46 Mean :356.67
## 3rd Qu.:24.000 3rd Qu.:666.0 3rd Qu.:20.20 3rd Qu.:396.23
## Max. :24.000 Max. :711.0 Max. :22.00 Max. :396.90
## lstat medv
## Min. : 1.73 Min. : 5.00
## 1st Qu.: 6.95 1st Qu.:17.02
## Median :11.36 Median :21.20
## Mean :12.65 Mean :22.53
## 3rd Qu.:16.95 3rd Qu.:25.00
## Max. :37.97 Max. :50.00Здесь мы можем видеть, что переменные «crim» и «black» принимают широкий диапазон значений.
Переменные «crim», «zn», «rm» и «black» имеют большую разницу между их медианной и средним значением, что указывает на множество выбросов в соответствующих переменных.
par(mfrow = c(1, 4))
boxplot(Boston$crim, main='crim',col='Sky Blue')
boxplot(Boston$zn, main='zn',col='Sky Blue')
boxplot(Boston$rm, main='rm',col='Sky Blue')
boxplot(Boston$black, main='black',col='Sky Blue')
Как и было предсказано, мы видим много выбросов.
Посчитаем корреляцию между признаками.
corr_matrix <- cor(Boston)
corrplot.mixed(corr_matrix)
Мы видим крайне сильную корреляцию - 0,9 - между индексом доступности радиальных магистралей (rad) и полной ставкой налога на имущество (tax). Также мы можем видеть сильную отрицательную корреляцию между расстоянием до центров занятости (dis) и долей некоммерческого бизнеса (indus), концентрацией оксидов азота (nox), долей старых домов (age); и между более низким статусом населения (lstat) и средней стоимостью домов, занимаемых владельцами (medv). И сильная положительная корреляция между indus и nox, indus и tax, nox и age, rm и medv.
Также мы видим мультиколлинеарность некоторых предикторов, что может стать проблемой при построении модели, поскольку даже небольшие изменения в любом из предикторов могут вызвать большие изменения в моделях.
4. Построение модели, прогнозирующей стоимость дома
Посмотрим распределение переменной medv.
qplot(x=medv,data=Boston,geom='histogram')## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
На данной гистограмме распределения средней цены дома мы видим, что данные немного смещены влево. Однако, линейную регрессию нельзя строить на ненормально распределенных данных. Чтобы исправить распределние прологарифмируем параметр medv.
qplot(x=log(medv),data=Boston,geom='histogram')## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
После логарифмирования у нас остался небольшой сдвиг вправо, но распределение выглядит более “нормальным”.
Задание 1. Построение полной модели.
Стандартизируем предикторы:
# сделаем наличие реки фактором
Boston$chas <- as.factor(Boston$chas)
stand_Boston <- Boston %>% mutate_at(c('crim', 'zn', 'indus', 'nox', 'age', 'dis', 'tax', 'ptratio', 'black', 'lstat'), scale)model1_stand = lm(log(medv)~., data= stand_Boston)
summary(model1_stand)##
## Call:
## lm(formula = log(medv) ~ ., data = stand_Boston)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.73361 -0.09747 -0.01657 0.09629 0.86435
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.320438 0.104407 22.225 < 2e-16 ***
## crim -0.088351 0.011315 -7.808 3.52e-14 ***
## zn 0.027345 0.012815 2.134 0.033349 *
## indus 0.016923 0.016886 1.002 0.316755
## chas1 0.100888 0.034486 2.925 0.003598 **
## nox -0.090199 0.017717 -5.091 5.07e-07 ***
## rm 0.090833 0.016728 5.430 8.87e-08 ***
## age 0.005928 0.014883 0.398 0.690567
## dis -0.103364 0.016811 -6.149 1.62e-09 ***
## rad 0.014267 0.002656 5.373 1.20e-07 ***
## tax -0.105466 0.025368 -4.157 3.80e-05 ***
## ptratio -0.082856 0.011337 -7.309 1.10e-12 ***
## black 0.037757 0.009815 3.847 0.000135 ***
## lstat -0.207344 0.014496 -14.304 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1899 on 492 degrees of freedom
## Multiple R-squared: 0.7896, Adjusted R-squared: 0.7841
## F-statistic: 142.1 on 13 and 492 DF, p-value: < 2.2e-16par(mfrow=c(2,2))
plot(model1_stand)
Полная модель:
R-squared составляет 0.7841
F-statistic составляет 142.1
Полная модель объясняет 78% изменчивости. Протестируем, есть ли параметры, которые незначительно влияют на модель. Для этого воспользуемся критерием Акаике.
step = stepAIC(model1_stand,direction = 'both')## Start: AIC=-1667.19
## log(medv) ~ crim + zn + indus + chas + nox + rm + age + dis +
## rad + tax + ptratio + black + lstat
##
## Df Sum of Sq RSS AIC
## - age 1 0.0057 17.755 -1669.0
## - indus 1 0.0362 17.786 -1668.2
## <none> 17.749 -1667.2
## - zn 1 0.1643 17.914 -1664.5
## - chas 1 0.3088 18.058 -1660.5
## - black 1 0.5339 18.283 -1654.2
## - tax 1 0.6235 18.373 -1651.7
## - nox 1 0.9351 18.684 -1643.2
## - rad 1 1.0413 18.791 -1640.3
## - rm 1 1.0637 18.813 -1639.7
## - dis 1 1.3639 19.113 -1631.7
## - ptratio 1 1.9270 19.676 -1617.0
## - crim 1 2.1995 19.949 -1610.1
## - lstat 1 7.3809 25.130 -1493.2
##
## Step: AIC=-1669.03
## log(medv) ~ crim + zn + indus + chas + nox + rm + dis + rad +
## tax + ptratio + black + lstat
##
## Df Sum of Sq RSS AIC
## - indus 1 0.0363 17.791 -1670.0
## <none> 17.755 -1669.0
## + age 1 0.0057 17.749 -1667.2
## - zn 1 0.1593 17.914 -1666.5
## - chas 1 0.3138 18.069 -1662.2
## - black 1 0.5431 18.298 -1655.8
## - tax 1 0.6205 18.376 -1653.7
## - nox 1 0.9645 18.720 -1644.3
## - rad 1 1.0356 18.791 -1642.3
## - rm 1 1.1452 18.900 -1639.4
## - dis 1 1.5471 19.302 -1628.8
## - ptratio 1 1.9224 19.677 -1619.0
## - crim 1 2.1988 19.954 -1612.0
## - lstat 1 8.1949 25.950 -1479.0
##
## Step: AIC=-1670
## log(medv) ~ crim + zn + chas + nox + rm + dis + rad + tax + ptratio +
## black + lstat
##
## Df Sum of Sq RSS AIC
## <none> 17.791 -1670.0
## + indus 1 0.0363 17.755 -1669.0
## + age 1 0.0058 17.786 -1668.2
## - zn 1 0.1451 17.936 -1667.9
## - chas 1 0.3399 18.131 -1662.4
## - black 1 0.5344 18.326 -1657.0
## - tax 1 0.6139 18.405 -1654.8
## - nox 1 0.9350 18.726 -1646.1
## - rad 1 1.0088 18.800 -1644.1
## - rm 1 1.1171 18.909 -1641.2
## - dis 1 1.7385 19.530 -1624.8
## - ptratio 1 1.8862 19.678 -1621.0
## - crim 1 2.2229 20.014 -1612.4
## - lstat 1 8.1604 25.952 -1481.0step$anova## Stepwise Model Path
## Analysis of Deviance Table
##
## Initial Model:
## log(medv) ~ crim + zn + indus + chas + nox + rm + age + dis +
## rad + tax + ptratio + black + lstat
##
## Final Model:
## log(medv) ~ crim + zn + chas + nox + rm + dis + rad + tax + ptratio +
## black + lstat
##
##
## Step Df Deviance Resid. Df Resid. Dev AIC
## 1 492 17.74938 -1667.194
## 2 - age 1 0.005723781 493 17.75510 -1669.031
## 3 - indus 1 0.036264380 494 17.79137 -1669.999step_model = lm(log(medv) ~ crim + zn + chas + nox + rm + dis + rad + tax + ptratio + black + lstat,data=stand_Boston)
summary(step_model)##
## Call:
## lm(formula = log(medv) ~ crim + zn + chas + nox + rm + dis +
## rad + tax + ptratio + black + lstat, data = stand_Boston)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.73400 -0.09460 -0.01771 0.09782 0.86290
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.328996 0.101070 23.043 < 2e-16 ***
## crim -0.088757 0.011298 -7.856 2.49e-14 ***
## zn 0.025361 0.012637 2.007 0.045308 *
## chas1 0.105148 0.034228 3.072 0.002244 **
## nox -0.083634 0.016414 -5.095 4.97e-07 ***
## rm 0.090673 0.016281 5.569 4.20e-08 ***
## dis -0.108878 0.015671 -6.948 1.18e-11 ***
## rad 0.013446 0.002540 5.293 1.82e-07 ***
## tax -0.094023 0.022774 -4.129 4.28e-05 ***
## ptratio -0.081025 0.011196 -7.237 1.77e-12 ***
## black 0.037679 0.009781 3.852 0.000133 ***
## lstat -0.204262 0.013570 -15.053 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1898 on 494 degrees of freedom
## Multiple R-squared: 0.7891, Adjusted R-squared: 0.7844
## F-statistic: 168.1 on 11 and 494 DF, p-value: < 2.2e-16Новая модель:
R-squared составляет 0.7864
F-statistic составляет 104.3
С помощью функции stepAIC мы последовательно рассчитали AIC для некоторых моделей и удалили предикторы, которые незначительно влияют на модель, что позволило несколько увеличить R-squared.
Задание 2. Проведем диагностику полной модели.
- Проверим на мультиколлениарность, которую мы заметили ранее.
library(car)
vif(step_model)## crim zn chas nox rm dis rad tax
## 1.789704 2.239229 1.059819 3.778011 1.834806 3.443420 6.861126 7.272386
## ptratio black lstat
## 1.757681 1.341559 2.581984Для дальнейшей работы нужно удалить предикторы с VIF > 5 (rad, tax).
- Проверим остатки:
plot(step_model$residuals)
abline(h=0,col='red')
Мы видим что остатки распределены ненормально.
- Найдем выбросы:
outlierTest(step_model)## rstudent unadjusted p-value Bonferroni p
## 413 4.775740 2.3647e-06 0.0011965
## 372 4.340845 1.7232e-05 0.0087195
## 373 4.007183 7.0941e-05 0.0358960
## 402 -3.952037 8.8814e-05 0.0449400Мы видим 4 выброса, влияющие на нашу модель.
Наибольший вклад вносит переменная lstat - -0.207344.
Задание 3. Постройте график предсказаний стоимости от переменной, которая обладает
наибольшим по модулю коэффициентом.
mod_hi <- lm(medv ~ as.numeric(lstat), data = stand_Boston)
test$predicted.medv <- predict(mod_hi,test)
pl1 <-test %>%
ggplot(aes(medv,predicted.medv)) +
geom_point(alpha=0.5) +
stat_smooth(aes(colour='black')) +
xlab('Значения medv') +
ylab('Предсказанные значения medv') +
theme_bw()
ggplotly(pl1)## `geom_smooth()` using method = 'loess' and formula 'y ~ x'Дополнительное задание.
Продолжим преобразовывать модель на НЕстандартизованных данных.
model1_stand = lm(log(medv)~., data= stand_Boston)
summary(model1_stand)##
## Call:
## lm(formula = log(medv) ~ ., data = stand_Boston)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.73361 -0.09747 -0.01657 0.09629 0.86435
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.320438 0.104407 22.225 < 2e-16 ***
## crim -0.088351 0.011315 -7.808 3.52e-14 ***
## zn 0.027345 0.012815 2.134 0.033349 *
## indus 0.016923 0.016886 1.002 0.316755
## chas1 0.100888 0.034486 2.925 0.003598 **
## nox -0.090199 0.017717 -5.091 5.07e-07 ***
## rm 0.090833 0.016728 5.430 8.87e-08 ***
## age 0.005928 0.014883 0.398 0.690567
## dis -0.103364 0.016811 -6.149 1.62e-09 ***
## rad 0.014267 0.002656 5.373 1.20e-07 ***
## tax -0.105466 0.025368 -4.157 3.80e-05 ***
## ptratio -0.082856 0.011337 -7.309 1.10e-12 ***
## black 0.037757 0.009815 3.847 0.000135 ***
## lstat -0.207344 0.014496 -14.304 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1899 on 492 degrees of freedom
## Multiple R-squared: 0.7896, Adjusted R-squared: 0.7841
## F-statistic: 142.1 on 13 and 492 DF, p-value: < 2.2e-16par(mfrow=c(2,2))
plot(model1_stand)
Новая модель:
R-squared составляет 0.7841
F-statistic составляет 142.1
step = stepAIC(model1_stand,direction = 'both')## Start: AIC=-1667.19
## log(medv) ~ crim + zn + indus + chas + nox + rm + age + dis +
## rad + tax + ptratio + black + lstat
##
## Df Sum of Sq RSS AIC
## - age 1 0.0057 17.755 -1669.0
## - indus 1 0.0362 17.786 -1668.2
## <none> 17.749 -1667.2
## - zn 1 0.1643 17.914 -1664.5
## - chas 1 0.3088 18.058 -1660.5
## - black 1 0.5339 18.283 -1654.2
## - tax 1 0.6235 18.373 -1651.7
## - nox 1 0.9351 18.684 -1643.2
## - rad 1 1.0413 18.791 -1640.3
## - rm 1 1.0637 18.813 -1639.7
## - dis 1 1.3639 19.113 -1631.7
## - ptratio 1 1.9270 19.676 -1617.0
## - crim 1 2.1995 19.949 -1610.1
## - lstat 1 7.3809 25.130 -1493.2
##
## Step: AIC=-1669.03
## log(medv) ~ crim + zn + indus + chas + nox + rm + dis + rad +
## tax + ptratio + black + lstat
##
## Df Sum of Sq RSS AIC
## - indus 1 0.0363 17.791 -1670.0
## <none> 17.755 -1669.0
## + age 1 0.0057 17.749 -1667.2
## - zn 1 0.1593 17.914 -1666.5
## - chas 1 0.3138 18.069 -1662.2
## - black 1 0.5431 18.298 -1655.8
## - tax 1 0.6205 18.376 -1653.7
## - nox 1 0.9645 18.720 -1644.3
## - rad 1 1.0356 18.791 -1642.3
## - rm 1 1.1452 18.900 -1639.4
## - dis 1 1.5471 19.302 -1628.8
## - ptratio 1 1.9224 19.677 -1619.0
## - crim 1 2.1988 19.954 -1612.0
## - lstat 1 8.1949 25.950 -1479.0
##
## Step: AIC=-1670
## log(medv) ~ crim + zn + chas + nox + rm + dis + rad + tax + ptratio +
## black + lstat
##
## Df Sum of Sq RSS AIC
## <none> 17.791 -1670.0
## + indus 1 0.0363 17.755 -1669.0
## + age 1 0.0058 17.786 -1668.2
## - zn 1 0.1451 17.936 -1667.9
## - chas 1 0.3399 18.131 -1662.4
## - black 1 0.5344 18.326 -1657.0
## - tax 1 0.6139 18.405 -1654.8
## - nox 1 0.9350 18.726 -1646.1
## - rad 1 1.0088 18.800 -1644.1
## - rm 1 1.1171 18.909 -1641.2
## - dis 1 1.7385 19.530 -1624.8
## - ptratio 1 1.8862 19.678 -1621.0
## - crim 1 2.2229 20.014 -1612.4
## - lstat 1 8.1604 25.952 -1481.0step$anova## Stepwise Model Path
## Analysis of Deviance Table
##
## Initial Model:
## log(medv) ~ crim + zn + indus + chas + nox + rm + age + dis +
## rad + tax + ptratio + black + lstat
##
## Final Model:
## log(medv) ~ crim + zn + chas + nox + rm + dis + rad + tax + ptratio +
## black + lstat
##
##
## Step Df Deviance Resid. Df Resid. Dev AIC
## 1 492 17.74938 -1667.194
## 2 - age 1 0.005723781 493 17.75510 -1669.031
## 3 - indus 1 0.036264380 494 17.79137 -1669.999step_model = lm(log(medv) ~ crim + zn + chas + nox + rm + dis + rad + tax + ptratio + black + lstat,data=stand_Boston)
summary(step_model)##
## Call:
## lm(formula = log(medv) ~ crim + zn + chas + nox + rm + dis +
## rad + tax + ptratio + black + lstat, data = stand_Boston)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.73400 -0.09460 -0.01771 0.09782 0.86290
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.328996 0.101070 23.043 < 2e-16 ***
## crim -0.088757 0.011298 -7.856 2.49e-14 ***
## zn 0.025361 0.012637 2.007 0.045308 *
## chas1 0.105148 0.034228 3.072 0.002244 **
## nox -0.083634 0.016414 -5.095 4.97e-07 ***
## rm 0.090673 0.016281 5.569 4.20e-08 ***
## dis -0.108878 0.015671 -6.948 1.18e-11 ***
## rad 0.013446 0.002540 5.293 1.82e-07 ***
## tax -0.094023 0.022774 -4.129 4.28e-05 ***
## ptratio -0.081025 0.011196 -7.237 1.77e-12 ***
## black 0.037679 0.009781 3.852 0.000133 ***
## lstat -0.204262 0.013570 -15.053 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1898 on 494 degrees of freedom
## Multiple R-squared: 0.7891, Adjusted R-squared: 0.7844
## F-statistic: 168.1 on 11 and 494 DF, p-value: < 2.2e-16Новая модель:
R-squared составляет 0.7844
F-statistic составляет 168.1
library(car)
vif(step_model)## crim zn chas nox rm dis rad tax
## 1.789704 2.239229 1.059819 3.778011 1.834806 3.443420 6.861126 7.272386
## ptratio black lstat
## 1.757681 1.341559 2.581984outlierTest(step_model)## rstudent unadjusted p-value Bonferroni p
## 413 4.775740 2.3647e-06 0.0011965
## 372 4.340845 1.7232e-05 0.0087195
## 373 4.007183 7.0941e-05 0.0358960
## 402 -3.952037 8.8814e-05 0.0449400У нас есть 4 выброса, которые мы удалим.
Boston_1 = Boston[-c(413,372, 373, 402),,drop=T]# drop the points
row.names(Boston_1)=1:nrow(Boston_1)
model2 = lm(log(medv)~crim+zn+chas+nox+dis+ptratio+black+lstat,data=Boston)
summary(model2)##
## Call:
## lm(formula = log(medv) ~ crim + zn + chas + nox + dis + ptratio +
## black + lstat, data = Boston)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.70805 -0.12147 -0.01373 0.10405 0.82894
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.7010244 0.1421502 33.071 < 2e-16 ***
## crim -0.0080105 0.0012719 -6.298 6.65e-10 ***
## zn 0.0013826 0.0005565 2.484 0.013311 *
## chas1 0.1302984 0.0363485 3.585 0.000371 ***
## nox -0.7322562 0.1352610 -5.414 9.63e-08 ***
## dis -0.0578514 0.0077684 -7.447 4.25e-13 ***
## ptratio -0.0376285 0.0048687 -7.729 6.06e-14 ***
## black 0.0002958 0.0001120 2.642 0.008505 **
## lstat -0.0353753 0.0017219 -20.545 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2025 on 497 degrees of freedom
## Multiple R-squared: 0.7584, Adjusted R-squared: 0.7545
## F-statistic: 195 on 8 and 497 DF, p-value: < 2.2e-16plot(model2)
Новая модель:
- R-squared составляет 0.7545
- F-statistic составляет 195
Не смотря на то, что мы устранили мульитколлинеанрость и незначимые предикторы, предсказания модели стали немного хуже (75%).
vif(model2)## crim zn chas nox dis ptratio black lstat
## 1.473811 2.074443 1.049493 3.024847 3.294717 1.367971 1.286382 1.861558Проведем диагностику модели
# Assumption - mean of the residuals is = 0
mean(model2$residuals)## [1] 4.05783e-18plot(model2$residuals)
abline(h=0,col='red')
Остатки распределены более меннее нормально.
Найдем выбросы.
outlierTest(model2)## rstudent unadjusted p-value Bonferroni p
## 413 4.270130 2.3419e-05 0.011850
## 372 4.052312 5.8859e-05 0.029782У нас есть два выброса, которые можно удалить.
Boston_1 = Boston[-c(413,372),,drop=T]# drop the points
row.names(Boston_1)=1:nrow(Boston_1)
model3 = lm(log(medv) ~ crim + zn + chas + nox + rm + dis +
ptratio + black + lstat,data=Boston_1)
par(mfrow=c(2,2))
plot(model3)
summary(model3)##
## Call:
## lm(formula = log(medv) ~ crim + zn + chas + nox + rm + dis +
## ptratio + black + lstat, data = Boston_1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.70878 -0.09807 -0.00993 0.08974 0.77161
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.6741127 0.1896496 19.373 < 2e-16 ***
## crim -0.0084826 0.0011748 -7.220 1.97e-12 ***
## zn 0.0007490 0.0005193 1.442 0.149874
## chas1 0.1199066 0.0335829 3.570 0.000391 ***
## nox -0.6147940 0.1258689 -4.884 1.40e-06 ***
## rm 0.1123588 0.0157923 7.115 3.96e-12 ***
## dis -0.0447458 0.0073176 -6.115 1.97e-09 ***
## ptratio -0.0333834 0.0045429 -7.348 8.36e-13 ***
## black 0.0004409 0.0001050 4.197 3.20e-05 ***
## lstat -0.0289046 0.0018797 -15.377 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1869 on 494 degrees of freedom
## Multiple R-squared: 0.7937, Adjusted R-squared: 0.7899
## F-statistic: 211.1 on 9 and 494 DF, p-value: < 2.2e-16Новая модель:
R-squared составляет 0.7899
F-statistic составляет 211.1
После удаления выбросов, мы видим, что значение R-squared увеличилось. Однако один из параметров перестал быть значимым.
model3_step=stepAIC(model3,Boston_1,direction = 'both')## Start: AIC=-1680.97
## log(medv) ~ crim + zn + chas + nox + rm + dis + ptratio + black +
## lstat
##
## Df Sum of Sq RSS AIC
## <none> 17.247 -1681.0
## - zn 1 0.0726 17.320 -1680.8
## - chas 1 0.4451 17.692 -1670.1
## - black 1 0.6151 17.862 -1665.3
## - nox 1 0.8329 18.080 -1659.2
## - dis 1 1.3054 18.552 -1646.2
## - rm 1 1.7673 19.014 -1633.8
## - crim 1 1.8202 19.067 -1632.4
## - ptratio 1 1.8853 19.132 -1630.7
## - lstat 1 8.2551 25.502 -1485.8model3_step$anov## Stepwise Model Path
## Analysis of Deviance Table
##
## Initial Model:
## log(medv) ~ crim + zn + chas + nox + rm + dis + ptratio + black +
## lstat
##
## Final Model:
## log(medv) ~ crim + zn + chas + nox + rm + dis + ptratio + black +
## lstat
##
##
## Step Df Deviance Resid. Df Resid. Dev AIC
## 1 494 17.24699 -1680.969summary(model3_step)##
## Call:
## lm(formula = log(medv) ~ crim + zn + chas + nox + rm + dis +
## ptratio + black + lstat, data = Boston_1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.70878 -0.09807 -0.00993 0.08974 0.77161
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.6741127 0.1896496 19.373 < 2e-16 ***
## crim -0.0084826 0.0011748 -7.220 1.97e-12 ***
## zn 0.0007490 0.0005193 1.442 0.149874
## chas1 0.1199066 0.0335829 3.570 0.000391 ***
## nox -0.6147940 0.1258689 -4.884 1.40e-06 ***
## rm 0.1123588 0.0157923 7.115 3.96e-12 ***
## dis -0.0447458 0.0073176 -6.115 1.97e-09 ***
## ptratio -0.0333834 0.0045429 -7.348 8.36e-13 ***
## black 0.0004409 0.0001050 4.197 3.20e-05 ***
## lstat -0.0289046 0.0018797 -15.377 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1869 on 494 degrees of freedom
## Multiple R-squared: 0.7937, Adjusted R-squared: 0.7899
## F-statistic: 211.1 on 9 and 494 DF, p-value: < 2.2e-16Новая модель:
R-squared составляет 0.7899
F-statistic составляет 211.1
par(mfrow=c(2,2))
plot(model3_step)
Мы видим, что нет разницы между model3 и model3_step, а также Р-квадрат имеет одинаковые значения. Теперь можем проанализировать выбросы модели.
outlierTest(model3)## rstudent unadjusted p-value Bonferroni p
## 369 4.303254 2.0307e-05 0.010235
## 372 4.266323 2.3833e-05 0.012012Мы видим 2 выбросы, которые мы можем выбросить и протестировать модель еще раз.
Boston_2=Boston_1[-c(371,400,373),,drop=T]
row.names(Boston_2)=1:nrow(Boston_2)
model4 = lm(log(medv) ~ crim + zn + chas + nox + rm + dis +
tax + ptratio + black + lstat,data=Boston_2)
summary(model4)##
## Call:
## lm(formula = log(medv) ~ crim + zn + chas + nox + rm + dis +
## tax + ptratio + black + lstat, data = Boston_2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.70835 -0.10150 -0.00960 0.08772 0.79154
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.624e+00 1.899e-01 19.086 < 2e-16 ***
## crim -7.880e-03 1.217e-03 -6.473 2.34e-10 ***
## zn 7.611e-04 5.218e-04 1.459 0.14530
## chas1 1.062e-01 3.335e-02 3.185 0.00154 **
## nox -5.614e-01 1.376e-01 -4.079 5.27e-05 ***
## rm 1.154e-01 1.550e-02 7.445 4.38e-13 ***
## dis -4.340e-02 7.194e-03 -6.033 3.18e-09 ***
## tax -6.026e-05 8.344e-05 -0.722 0.47055
## ptratio -3.255e-02 4.937e-03 -6.594 1.11e-10 ***
## black 4.396e-04 1.048e-04 4.195 3.23e-05 ***
## lstat -2.858e-02 1.870e-03 -15.288 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1831 on 490 degrees of freedom
## Multiple R-squared: 0.7971, Adjusted R-squared: 0.7929
## F-statistic: 192.4 on 10 and 490 DF, p-value: < 2.2e-16Новая модель:
R-squared составляет 0.7929
F-statistic составляет 192.4
par(mfrow=c(2,2))
plot(model4)
Мы видим, что мы снова увеличили мощность нашей модели, но потеряли значение по другой переменной. Мы можем запустить stepAIC и посмотреть, сможем ли мы его удалить.
outlierTest(model4)## rstudent unadjusted p-value Bonferroni p
## 369 4.527033 7.5217e-06 0.0037683
## 371 4.517479 7.8548e-06 0.0039352
## 372 3.982469 7.8571e-05 0.0393640
## 398 -3.953793 8.8281e-05 0.0442290Мы видим 4 выбросы, которые мы можем выбросить и протестировать модель еще раз.
Boston_3 =Boston_2[-c(369,371,372,398),,drop=T]
row.names(Boston_3)=1:nrow(Boston_3)
Boston_3$chas <- as.factor(Boston_3$chas)
model5 = lm(log(medv) ~ crim + zn + chas + nox + rm + dis +
tax + ptratio + black + lstat,data=Boston_3)
summary(model5)##
## Call:
## lm(formula = log(medv) ~ crim + zn + chas + nox + rm + dis +
## tax + ptratio + black + lstat, data = Boston_3)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.64945 -0.09613 -0.01079 0.08989 0.69041
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.421e+00 1.781e-01 19.207 < 2e-16 ***
## crim -7.826e-03 1.132e-03 -6.915 1.48e-11 ***
## zn 5.723e-04 4.854e-04 1.179 0.239025
## chas1 7.966e-02 3.144e-02 2.533 0.011610 *
## nox -5.005e-01 1.282e-01 -3.905 0.000108 ***
## rm 1.385e-01 1.468e-02 9.436 < 2e-16 ***
## dis -3.810e-02 6.719e-03 -5.670 2.45e-08 ***
## tax -1.143e-04 7.801e-05 -1.465 0.143473
## ptratio -3.225e-02 4.591e-03 -7.023 7.35e-12 ***
## black 4.324e-04 9.785e-05 4.418 1.23e-05 ***
## lstat -2.677e-02 1.789e-03 -14.961 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1701 on 486 degrees of freedom
## Multiple R-squared: 0.82, Adjusted R-squared: 0.8163
## F-statistic: 221.4 on 10 and 486 DF, p-value: < 2.2e-16Новая модель:
R-squared составляет 0.8163
F-statistic составляет 221.4
par(mfrow=c(2,2))
plot(model5)
Мы улучшили еще немного нашу модель и теперь она объясняет 82% изменчивости. Нет предела совершенству, меня устраивает модель, объясняющая 82% изменчивости. Наибольший вклад вносят параметры: количество комнат, количество оксида азота, отношение учеников к учителям. Так что я бы посоветовала клиенту, строить большие дома рядом со школой и с парком неподалеку.
Финальный анализ модели.
residuals <- data.frame('Residuals' = model5$residuals)
res_hist <- ggplot(residuals, aes(x=Residuals)) + geom_histogram(color='black', fill='skyblue') + ggtitle('Histogram of Residuals')
res_hist## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Остатки распределены нормально.
plot(model5, col='Sky Blue')
# пересчитаем тестовый датасет
set.seed(123)
split <- sample.split(Boston,SplitRatio = 0.75) #assigns booleans to a new coloumn based on the split ratio
train <- subset(Boston,split==TRUE)
test <- subset(Boston,split==FALSE)
test$predicted.medv <- predict(model5,test)
pl1 <-test %>%
ggplot(aes(medv,predicted.medv)) +
geom_point(alpha=0.5) +
stat_smooth(aes(colour='black')) +
xlab('Значения medv') +
ylab('Предсказанные значения medv') +
theme_bw()
ggplotly(pl1)## `geom_smooth()` using method = 'loess' and formula 'y ~ x'